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Abstract 

We examined individual heterogeneity in survival and recruitment of female 
Pacific black brant (Branta bemicla nigricans) using frailty models adapted to a 
capture-mark-recapture context. Our main objectives were (1) to quantify lev- 
els of heterogeneity and examine factors affecting heterogeneity, and (2) model 
the effects of individual heterogeneity on harvest dynamics through matrix 
models. We used 24 years of data on brant marked and recaptured at the Tuta- 
koke River colony, AK. Multievent models were fit as hidden Markov chain 
using program E-SURGE with an adequate overdispersion coefficient. Annual 
survival of individuals marked as goslings was heterogeneous among individuals 
and year specific with about 0.23 difference in survival between "high" (0.73)- 
and "low" (0.50) -quality individuals at average survival probability. Adult 
survival (0.85 ± 0.004) was homogeneous and higher than survival of both 
groups of juveniles. The annual recruitment probability was heterogeneous for 
brant >l-year-old; 0.56 (±0.21) and 0.31 (±0.03) for high- and low-quality 
individuals, respectively. Assuming equal clutch sizes for high- and low-quality 
individuals and that 80% of offspring were in the same quality class as the 
breeding female resulted in reproductive values about twice as high for high- 
quality individuals than low-quality individual for a given class of individuals 
producing differential contributions to population growth among groups. Dif- 
ferences in reproductive values greatly increased when we assumed high-quality 
individuals had larger clutch sizes. When we assumed that 50% of offspring 
were in the same quality class as their mothers and clutches were equal, 
differences in reproductive values between quality classes were greatly reduced 
or eliminated (breeders [BRs]). We considered several harvest scenarios using 
the assumption that 80% of offspring were in the same quality class as their 
mothers. The amount of compensation for harvest mortality declined as the 
proportion of high-quality individuals in the harvest increased, as differences in 
clutch sizes between groups decreased and as the proportion of BRs in the har- 
vest increased. Synthesis and applications. Harvest at the same proportional level 
of the overall population can result in variable responses in population growth 
rate when heterogeneity is present in a population, k was <1.0 under every sce- 
nario when harvest rates were >10%, and heterogeneity caused as much as +2% 
difference in growth rates at the highest levels of proportional harvest for low- 
quality individuals and the greatest differences in qualities between classes of 
individuals, a critical difference for a population with ). near 1.0 such as the 
brant. We observed less response in overall survival in the presence of hetero- 
geneity because we did not observe heterogeneity in the annual survival of BRs. 
This analysis provides a comprehensive view of overall compensation at the 
population level and also constitutes the first example of a survival-recruitment 
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model with heterogeneity. Individual heterogeneity should be more explicitly 
considered in harvest management of vertebrates. 



Introduction 

The harvest of wildlife and fisheries populations has been 
the subject of considerable debate for almost a century 
(Baranov 1918; Beverton and Holt 1957). Difficulty with 
monitoring fish stocks and unexpected changes in harvest 
encouraged the use of modeling (e.g., Getz and Haight 
1989) and the development of theory to guide harvest 
management, which was then applied to a number of 
North American waterfowl populations (Anderson and 
Burnham 1976). More recently, the issue of harvest was 
examined as a conservation problem (Reynolds et al. 
2001), and symmetrically incidental exploitation (e.g., dif- 
fuse mortality induced by human activities in an other- 
wise protected species) was considered as an exploitation 
problem (Lebreton 2005). Whatever the context, a central 
question in harvest dynamics is that of compensation; 
does harvest, as a specific cause of mortality, add its 
effects totally independently of natural mortality or is the 
effect partially or totally compensated? 

Harvest management of most waterfowl, particularly in 
North America, is currently guided by the assumption 
that harvest mortality is at least partially additive (John- 
son et al. 1993; Conn and Kendall 2004), and if compen- 
sation occurs, it is primarily through density dependence 
in survival probability and reproduction (Anderson and 
Burnham 1976; Nichols et al. 1995). Harvest mortality 
may be compensated through density-dependent increases 
in survival or reproduction postharvest, such that harvest 
mortality may have no effect on overall survival or 
growth rate of the population. The evidence for additive 
or compensatory harvest mortality is mixed (Nichols 
et al. 1984; Rexstad 1992; Smith and Reynolds 1992; 
Gauthier et al. 2001; Williams et al. 2002; Sedinger et al. 
2007; Sedinger and Herzog 2012) and is least compen- 
sated in populations with higher inherent survival proba- 
bility. Evidence for density-dependent regulation of 
survival and reproduction in waterfowl is inconsistent 
(Johnson et al. 1992; Anderson et al. 1997; Sedinger et al. 
1998; Viljugrein et al. 2005) and may be too weak for 
compensation to occur (Lebreton 2005). Moreover, 
detecting density dependence (Lebreton 2009) or a nega- 
tive correlation between natural mortality and harvest 
(Schaub and Lebreton 2004) is the subject of notable sta- 
tistical difficulties (Otis and White 2004). As a conse- 
quence, some retrospective analyses may overestimate the 
prevalence of density dependence (Shenk et al. 1998; 
Lebreton 2005). Furthermore, distinguishing between the 
effects of harvest and density on abundance is difficult 
because harvest regulations are typically liberal when 



abundance is high and conservative when abundance is 
low (Smith and Reynolds 1992; Sedinger and Rexstad 
1994; Sedinger and Herzog 2012). 

Because of the continuing debate about the role of den- 
sity dependence in compensation, alternative functional 
forms of the relationship between harvest and population 
change may be required to fully represent these dynamics 
(Runge and Johnson 2002; Conn and Kendall 2004). We 
examined how heterogeneity or individual variation in 
survival and recruitment of Pacific black brant (Branta 
bernicla nigricans; hereafter brant) may provide an alter- 
nate explanation for the relationship between harvest and 
population dynamics. Others (Johnson et al. 1984, 1988; 
Lebreton 2005) previously demonstrated through statisti- 
cal and population models that individual heterogeneity 
in survival and reproduction can lead to compensation 
even in the absence of density dependence, and we fur- 
ther these findings using data on a specific population 
and new modeling tools. 

The consequences of harvest are linked to the expected 
contributions of the harvested individuals to future popu- 
lation growth, which are measured by reproductive value 
(MacArthur 1960; Kokko 2001). Harvest of an individual 
with high reproductive value has more effect on popula- 
tion dynamics than harvest of an individual of low repro- 
ductive value (Brooks and Lebreton 2001; Kokko 2001; 
Hauser et al. 2006). Disproportionate harvest risk of indi- 
viduals in poor physiological condition (Greenwood et al. 
1986; Dufour et al. 1993) is consistent with the hypothe- 
sis that heterogeneity in reproductive value provides some 
compensation for harvest mortality. Heterogeneity in 
reproductive value of individuals has been clearly linked 
to a number of characteristics (e.g., age, gender); how- 
ever, unexplained sources of heterogeneity may have addi- 
tional effects on reproductive value and harvest dynamics, 
and these sources of heterogeneity have created modeling 
challenges (Vaupel and Yashin 1985; Link et al. 2002). 

In human health studies, individual heterogeneity is 
considered as a random variable with a continuous distri- 
bution in so-called "frailty" survival models (Vaupel 
1990), and these models have also been successfully used 
in studies of free-ranging vertebrate (e.g., Cam et al. 
2002). An alternative is to consider that demographic 
parameters are distributed according to a discrete distri- 
bution; the population is a mixture of several types of 
individuals differing in parameters such as survival proba- 
bility. Pledger and Schwarz (2002) reviewed capture- 
recapture survival models with heterogeneity either 
continuous or based on finite mixtures, and they con- 
cluded that two-level mixtures (i.e., considering the 
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population is composed of two types of individuals) often 
provided an adequate representation of individual hetero- 
geneity (but see Dorazio and Royle 2003). We explored a 
novel application of finite mixture models to multiply 
demographic parameters. 

Rich, longitudinal data sets are still needed to model 
individual heterogeneity. The long-term study of brant in 
western Alaska (e.g., Sedinger et al. 2008; Nicolai et al. 
2012) presented an excellent opportunity for us to exam- 
ine new methods for simultaneously modeling heteroge- 
neity in capture, recruitment, and survival probability and 
explore how this variation may change through time. We 
use estimated levels of heterogeneity to explore implica- 
tions of this heterogeneity for understanding effects of 
historical harvest rates for brant on population dynamics. 
Our main objectives were (1) to quantify levels of hetero- 
geneity and examine factors affecting heterogeneity and 
(2) model the effects of individual heterogeneity on har- 
vest dynamics through matrix models. We predicted that 
individual heterogeneity in survival and recruitment 
would exist in the brant population, and that heterogeneity 
would compensate for some harvest. 

Materials and Methods 

Study system and data collection 

Data were collected at the Tutakoke River brant colony 
(61°159'N, 165°37W) on the Yukon-Kuskokwim Delta, 
AK, which historically represented about 20% of the breed- 
ing population (Sedinger et al. 1993). Brant are long- 
distance migrants that winter along the Pacific coast of 
North America from Alaska to the Pacific coast of Mexico 
(Reed et al. 1998). Brant from the Tutakoke River colony 
nest in coastal tundra within 2 km of the Bering Sea coast. 
We captured individual brant by herding them into corral 
traps during the adult remigial molt in mid-late July each 
summer (Sedinger et al. 1997). We assigned individuals to 
one of three age classes upon initial capture based on plum- 
age characteristics: goslings (HY, approximately 1 month 
old); 1-year-olds (SY,13 months old) and adults 
(ASY, >25 months old). We placed U.S. Geological Survey 
steel bands and uniquely engraved plastic bands on captured 
individuals; plastic bands could be read at a distance of 
>100 m with spotting scopes, which facilitated detection in 
subsequent years. We determined gender by cloacal exami- 
nation. Goslings in nests associated with marked adults 
received uniquely numbered webtags while still in the nest 
(Sedinger and Flint 1991), which allowed us to determine 
ages of goslings (±1 day) when captured in banding drives. 
We weighed and measured all web-tagged goslings when 
captured during banding drives. We calculated a cohort 
mean mass adjusted for gender and age (days). 



We restricted the analysis to females because (1) most 
males disperse to other breeding locations (Lindberg et al. 
1998), reducing the number of known-age individuals 
available for study; and (2) brant maintain long-term pair 
bonds, so fates of males and females were not indepen- 
dent. We encountered marked individuals during the 
aforementioned banding drives and during nesting. 
However, few 1 -year-olds were observed and none have 
ever been confirmed as breeders (BRs), encounters of 
1 -year-olds were therefore removed and we modeled 
initial survival over the first 2 years of life. At nesting, we 
visited 50, randomly located 50-m radius plots every 
4 days during egg laying and at least on alternate days 
during the hatch period. We read plastic bands when 
females were flushed from their nests. Additionally, we 
attempted to check all nesting females on the colony for 
plastic bands by flushing them from their nests, beginning 
shortly after the end of egg laying. We detected about 
79% of breeding females that were present on the colony 
annually (Sedinger et al. 2001). 

Multievent models 

Analysis was based on female brant marked with individ- 
ually coded tarsal tags and legbands at the Tutakoke River 
colony, AK from 1987 to 2009 (Sedinger et al. 1998). We 
examined individual heterogeneity of female brant based 
on CMR recruitment models (see e.g., Lebreton et al. 
2003) incorporating individual heterogeneity in the vari- 
ous parameters of these models. The resulting model was 
a multievent model (Pradel 2005), which is a type of 
hidden Markov chain model (Choquet et al. 2009). We 
considered a suite of models both with and without indi- 
vidual heterogeneity and other sources of variation to 
estimate state (i.e., quality) of an individual and capture, 
apparent survival, and recruitment probability. We used 
an information-theoretic approach for model selection 
and to determine what sources of variation were sup- 
ported (Burnham and Anderson 2002). When included, 
individual heterogeneity was modeled using a two-level 
mixture model with two hidden groups. This structure 
results in four states: 1 = BR group A, 2 = prebreeder 
(PB) group A, 3 = BR group B, 4 = PB group B. A short- 
coming of this approach is that heterogeneity in capture 
probability is linked with heterogeneity in survival and 
recruitment, such that an individual in a specific survival 
or recruitment group has to be in the same group for 
detection probability. Separating these two types of heter- 
ogeneity would require four hidden groups or eight states. 
Such a model would be quite unstable and present severe 
identifiability problems and was not considered. The neg- 
ative sampling correlations between parameters, however, 
will tend to exaggerate heterogeneity in the demographic 
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parameters of interest, and we consider that issue in the 
Discussion. 

As we are only sampling at a breeding colony, individ- 
uals are observed as a BR (only after second year) or PB 
(in their hatch year) at initial capture and as a BR in any 
subsequent recapture. However, we do not know to which 
heterogeneity state (A or B) an individual belongs because 
this is a hidden state. Therefore, the proportion of indi- 
viduals in heterogeneity groups A and B have to be esti- 
mated within each breeding category. These initial-state 
probabilities come as extra parameters in multievent 
models when compared to usual multistate CMR models, 
in which the state of an individual at time of marking is 
known. A model with year-specific variation in initial 
state and 25 occasions would require 75 parameters 
because the 4th state is estimated as 1— the sum of the 
probability of the other three states. 

Recapture probability (pi) was the probability that an 
individual alive and in the Tutakoke River population in 
year i was reencountered in that year. We considered 
models with effects of year and heterogeneity on p, with 
an additive and interactive relationship between these fac- 
tors. With p and other parameters, we considered models 
with linear trends (logit scale) in heterogeneity to exam- 
ine whether heterogeneity changed through years because 
the size of the brant colony was declining during the 
study period, which may have changed the amount of 
heterogeneity through time. Apparent survival probability 
(([>;), the probability that a brant alive and associated with 
the Tutakoke colony at year survives and does not 
permanently emigrate between year i and i + 1, was mod- 
eled as a function of age, year, mass at capture (i.e., mean 
year cohort mass adjusted for age; Sedinger and Nicolai 
2011), and heterogeneity. Juvenile and adult survival 
probabilities were modeled independently. Juvenile sur- 
vival spans the 2 years from fledging to potential age of 
first breeding, since second year encountered were 
removed. We estimated survival for the first and second 
year of life assuming survival was equal for these two 
annual periods. Recruitment probability (aj), the 
probability that animal of age j starts to breed at that age 
(Pradel and Lebreton 1999), was modeled as a function 
of age from 2-years old up to 6-years old with and 
without the effects of heterogeneity. 

We assessed fit of the model to the data using tests 
conducted in program U-CARE following Crespin et al. 
(2006) for recruitment models and Fletcher et al. (2012) 
for models with heterogeneity. By definition, heterogene- 
ity in capture probability corresponds to a mix of individ- 
uals with low and high capturability. Compared with 
homogeneity, the data set will present an excess of runs 
of 0 and of runs of 1. This situation results in the simul- 
taneous presence of various degrees of transience and trap 
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happiness, respectively, detected by goodness-of-fit (GOF) 
components TEST3.SR and TEST2.CT, as shown by 
Fletcher et al. (2012). We expect all other components to 
be somewhat sensitive to heterogeneity in all parameters. 

All models were run with random initial values 
repeated 10 times to protect against local minima. In all 
cases, the minimal value of the deviance was that 
obtained for the largest number of repeats, so we are con- 
fident that the deviances used in QAIC calculation corre- 
sponded to the global minimum. For a few of the most 
complex models, the deviance was obtained after 200 iter- 
ations and full convergence was not yet attained, implying 
a slight overestimation not bearing any consequence on 
QAIC-based model selection. This was always due to the 
slow optimization of boundary estimates (i.e., a probabil- 
ity estimated to be 1.0, has to converge in logit to "infin- 
ity," represented by the value 15 in E-SURGE), and this 
convergence requires many iterations, without substantial 
changes to the deviance. Our checks with a few of these 
models showed that the difference in deviance was always 
<0.100, far from the difference in QAIC with the pre- 
ferred models. 

Matrix models 

Using estimates from the best multievent models and 
published estimates of fecundity (Sedinger et al. 1998; 
Nicolai and Sedinger 2012), we parameterized a matrix 
model with two population segments (A and B) in a gos- 
ling, prebreeding, or breeding state (six states total) using 
a postbirth pulse structure (Caswell 2001). Matrix models 
were run in ULM (Legendre and Clobert 1995) and 
MATLAB® (Mathworks Inc., Natick, MA). We analyzed 
the matrix to determine the reproductive value of various 
states and changes in growth rate (A) under harvest sce- 
narios that included different proportions of brant from 
the six states. We calculated reproductive values in the 
presence of heterogeneity with different levels of inheri- 
tance of the states where differences in "heritability" rep- 
resented individuals producing different proportion of 
goslings in their quality state. We use the term inheri- 
tance rather than heritability, to acknowledge that both 
genetic and environmental (e.g., timing of reproduction, 
brood-rearing area) factors may affect quality of offspring. 
We also considered how reproductive values and X chan- 
ged with differences in clutch size among quality classes. 

Under the assumptions that females produced 80% off- 
spring of the same class, and in addition to heterogeneity 
in survival and recruitment, we evaluated the effects of 
harvest with different proportions of quality classes in the 
harvest, different clutch sizes between classes, and differ- 
ent proportions of goslings, PBs, or BRs in the harvest. 
Scenarios for harvest rates were developed based on 
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estimates from Sedinger et al. (2007). Harvest effects were 
evaluated by comparing resulting X to the per capita 
growth rate expected under a uniform harvest proportion 
(h) irrespective of reproductive value: 

X(U) = a(0)(1 - h); 

where i(0) is the growth rate in the absence of harvest 
and h is the overall proportion of the population in the 
harvest. When reproductive values and proportion of 
quality classes in the harvest are considered an approxi- 
mation of X is obtained based on sensitivity analysis as: 

X{H) « A(0)(1 - bh); 

where b is the ratio of harvest proportion weighted by 
reproductive value of a quality class (h n ) and harvest 
proportion irrespective of reproductive value (h). Harvest 
proportion irrespective of reproductive value is: 

n 

E w i h i 
E w i 

and harvest proportion weighted by reproductive value is: 

n 

E nnh 

n„ = —„ 

E w i v i 

i=1 

where w ; is proportion of population in age/quality class 
at stable structure and v,- is the reproductive value for the 
;th age/quality class. The effects of harvest on X when 
considering reproductive value will be smaller than the 
effects of uniform harvest irrespective of reproductive 
value when b < 1.0 (i.e., compensation occurs because 
harvest is proportionally higher in the population seg- 
ments with lower reproductive value) and b can therefore 
be interpreted as the relative strength of harvest effects 
with stronger effects and less compensation as b 
approaches 1.0 (i.e., if b = 0 then harvest would hypo- 
thetically have no effect and if b = 1.0 no compensation 
would occur). We also examined the effects of harvest on 
overall population-level survival for the harvest scenario 
with the maximum effects of heterogeneity because 
survival is the parameter often examined to evaluate 
compensation. 

Results 

Multievent models 

We included 31,167 individual female brant in our analy- 
sis and identified 2870 different capture histories. GOF 
tests provide strong evidence of heterogeneity. Following 



the results in Fletcher et al. (2012, pp. 208-209), 
components of TEST3.SR (z = 55.1) and TEST2.CT 
(z = —14.61), distributed as N(0, 1) under full homo- 
geneity, indicated strong heterogeneity in capture proba- 
bility, with, as expected a higher sensitivity of TEST3.SR. 
Other tests (nondirectional components of TEST3.SR and 
TEST2.CT + TEST3.SM and TEST2.CL, %j 9l = 919.5, 
P = 0.000) indicated heterogeneity in survival probability. 
If all test components were considered, the variance infla- 
tion factor (c) would equal 7.23 and 3.16 without the 
directional components directly linked to recapture heter- 
ogeneity. This latter value was used as an overdispersion 
factor in QAIC-based model selection because we consid- 
ered only models that included heterogeneity in P. 
Because of the sensitivity of the GOF components used to 
estimate c to heterogeneity in all parameters, we suspect 
that this c estimate is still biased high (White 2002); how- 
ever, we are unsure what lower value this estimate should 
be adjusted to and prefer instead to use a conservative 
approach to model selection (see also Lebreton et al. 
2003). 

A general model with full-year effects and heterogeneity 
had 338 identifiable parameters, QAIC = 41903.7, and 
served as a benchmark for model selection. Removal of 
heterogeneity from proportion of individuals in initial 
states, recruitment, and survival parameters, but leaving 
heterogeneity in recapture probability resulted in a model 
with the highest QAIC value (42557.301, for 123 identifi- 
able parameters) of the 65 models considered, thus con- 
firming the presence of strong demographic heterogeneity 
as observed from the goodness-of-fit tests, and not just 
heterogeneity in recapture probability. 

We found strong support for removing year effects in 
recruitment and the interaction between year and age 
effects on survival probability (AQAIC = —275.8 for 
model without these effects). Year trending in survival 
probability was also not supported (AQAIC = +29.2 for 
model with trending). A model with a different, year- 
specific variation in HY and ASY survival, and no hetero- 
geneity in ASY survival were preferable (AQAIC = —20.7) 
to a model with additive year effects for HY and ASY 
survival and heterogeneity across all age classes. 

Age-specific variation in recruitment probability was 
not strongly supported as AQAIC = 0.5 between a model 
with no age-specific transition to breeding (105 para- 
meters) and a model (107 parameters) with age-specific 
recruitment for age 2 and all other age classes (3+). 
Furthermore, the base model with age-specific recruit- 
ment for ages 2-5 had a AQAIC = +280.2 compared with 
the model without age-specific recruitment. Assuming a 
recruitment probability constant over age induces a 
cumulative recruitment curve slowly declining over years 
and ensures identifiability of the model. Models with full 
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recruitment completed at age 5, 6, or 7 did not do better 
than these models with constant recruitment. 

Models with additive effects of year and heterogeneity 
on capture probability generally had lower QAIC values 
than models with interactions between these factors as the 
top 17 models had additive effects. However, this result is 
somewhat in contrast to variability in components of 
TEST3.SR, which suggests a change in heterogeneity over 
time, the variability in the proportion of low-capturability 
individuals translating in an excess of individuals never 
seen again after first detection as a BR. For this same 
reason, we suspect that TEST3.SR components, which 
compare the proportion of individuals seen again or not 
seen again between newly marked and already marked 
individuals, may be sensitive to heterogeneity in initial- 
state probabilities rather than only to heterogeneity in 
recapture probability. 

The best approximating model included year-specific 
variation in probability that an individual was in group A 
or B at initial capture (Fig. 1) and was favored over the 
second (AQAIC = 3.7) and third best (AQAIC = 10.4) 
models, which included an additive relationship between 
proportion BR and PB and year or mass effects on HY 
survival, respectively. Therefore, we used estimates only 
from the top model. In this model, HY survival was 
modeled with an additive relationship between year and 
heterogeneity (Fig. 2) with differences in survival between 
high (A) and low (B) quality groups equal to 1.38 
(SE = 0.84) on a logit scale or about 0.23 at average (0.50 
low quality and 0.73 high quality) survival probability. 
ASY survival was constant and equaled 0.85 (SE = 0.004). 
Recruitment probability was constant across ages, but the 
rate of recruitment was higher for group A (0.54, 



SE = 0.214) than group B (0.31, SE = 0.028). Capture 
probabilities varied across years and included heterogene- 
ity in an additive relationship. Capture probabilities were 
lower for group A than B and may indicate a sampling 
covariance that exaggerates heterogeneity in other para- 
meters (see Discussion). 

Matrix models 

Assuming 80% of the goslings produced were in the same 
quality class as the breeding female and equal clutch size 
between quality classes resulted in a projected X = 1.058 
(Table 1). Under these conditions, reproductive values of 
high-quality brant were about twice as high as the com- 
parable states for low-quality brant. Differences in repro- 
ductive value increased considerably when we assumed 
the differences in quality also resulted in a difference in 
clutch size (5 vs. 3) and A increased to 1.090 (Table 1). 
When we assumed no inheritance of quality (i.e., 50% of 
offspring in class) and equal clutch size, differences in 
reproductive values between high- and low-quality 
goslings and PBs were less than under the other conditions 
and equal among classes of BRs because of identical ASY 
survival in this class of birds (Table 1). High-quality BRs 
had the highest reproductive values and low-quality goslings 
the lowest reproductive values under every scenario. 

The compensation index (b) increased (compensation 
decreased), and the relative effects of harvest therefore 
increased when (1) differences in clutch size between 
classes decreased, (2) the proportion of high-quality indi- 
viduals in the harvest increased, and (3) the proportion 
of BRs in the harvest increased. The compensation index 
was below 1.0 for all the scenarios we considered because 




oooooocncncncncncncnaiaicn 



oSSooooooS Figure 1. Year-specific probability that an 

^oooooooooo ..... , . . . . ■ ... i 

HNMNNNnNMNNN individual was in state A at initial capture. 

Year Error bars are ±1 standard error. 
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Table 1. Population growth rates and reproductive values of brant for the six states under different scenarios of inheritance and heterogeneity. 





Clutch size 


Clutch size 




Reproductive value 










Inheritance 


(high quality) 


(low quality) 


X 


GOS-high 


PB-high 


B-high 


GOS-low 


PB-low 


BR-low 


0.8 


4 


4 


1.058 


0.174 


0.251 


0.292 


0.046 


0.097 


0.140 


0.5 


4 


4 


1.044 


0.135 


0.192 


0.222 


0.074 


0.156 


0.222 


0.8 


5 


3 


1.090 


0.234 


0.346 


0.408 


0.002 


0.004 


0.006 



GOS, gosling; PB, prebreeders; BR, breeders; high, high-quality individuals; low, low-quality individuals. 



the ratio of low-quality individuals in the harvest <0.51 
(Table 2, Fig. 3). Under all scenarios and levels of com- 
pensation, 1 was below 1.0 when overall harvest rate was 
>10%. Under scenario 3 and at 10% overall harvest, 
compensation through heterogeneity resulted in 
a = 1.009; 0.028 higher (about 3% higher growth rate) 
compared with the situation when heterogeneity was not 
considered (X = 0.981). We observed almost no response 
in population survival under the harvest scenario 3 
because at stable structure, most of the population was 
composed of BRs, which exhibited no heterogeneity in 
survival (Fig. 4). 

Discussion 

Heterogeneity modeling 

Our analysis provided strong evidence that heterogeneity 
in survival and recruitment was present in the brant pop- 
ulation. To our knowledge, this is the first time that the 



levels of heterogeneity were quantified simultaneously for 
>1 parameter, and our approach may be useful for others 
examining heterogeneity. However, we acknowledge that 
additional advances would have improved our modeling. 
We did not simultaneously model heterogeneity in the 
full suite of population parameters (e.g., clutch size), and 
state space modeling (e.g., Borysiewicz et al. 2009) may 
be a useful tool for heterogeneity modeling when different 
types of data are available (e.g., longitudinal and survey 
data). Also, frailty models provided adequate structure for 
our questions, but modeling heterogeneity as a distribu- 
tion may be preferable for other questions (Dorazio and 
Royle 2003). Finally, we could not resolve the link between 
estimates of capture probability and quality classes 
with higher-quality individuals having lower-capture 
probabilities. Collectively, we do not think these issues 
had much effect on our estimates, except to possibly 
inflate levels of heterogeneity in survival and recruitment; 
however, we believe our inferences are relevant to the 
population of interest. 
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Table 2. Three scenarios of harvest for brant and resulting levels of compensation (b). 



Scenario 


Clutch size (high quality) 


Clutch size (low quality) 


Harvest ratio high/low 


Harvest ratio pre-breeder/breeder 


b 


1 


4 


4 


0.5 


2 


0.812 


2 


5 


3 


0.5 


2 


0.791 


3 


5 


3 


0.3 


3 


0.635 




0 06 



010 




0 06 



0 10 




0.02 0.04 0.06 0.08 
Proportion harvested 



010 



Figure 3. Population growth rate as a 
function of the proportion of the population 
harvested under scenario 1 (A), 2 (B), or 3 (C; 
Table 2). Thin dotted line = growth rate 
without considering heterogeneity (no 
compensation); thick solid line = growth rate 
with compensation for heterogeneity; thin 
solid line = growth rate with compensation for 
harvest parameter sensitivity. 



Sources of heterogeneity 

We observed yearly variation in the proportion of individ- 
uals in the initial quality states. We suspect this heteroge- 
neity is related to the quality of individuals breeding in a 



given year, the quality of the environment that goslings 
experienced prior to fledging, and genetic-environmental 
covariance (Sedinger and Chelgren 2007; Nicolai and 
Sedinger 2012). However, we are unclear why average mass 
of a cohort, which reflects quality of environment goslings 
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Figure 4. Changes in overall population 
survival (S) relative to changes in proportion of 
population harvested in a brant population 
without (dashed line) and with heterogeneity 
(solid line) as described under harvest scenario 
3 in Table 2. 



0.6 (- 
0.5 



0.02 0.04 0.06 0.08 0.1 0.12 
Proportion harvested 



0.14 0.16 



experience prior to fledging, did not receive additional 
support in our analysis. Perhaps, the strength of the three 
causes of heterogeneity identified above varies through 
time, and a single description of these sources (e.g., gosling 
mass) is an inadequate description (Sedinger and Chelgren 
2007). 

Heterogeneity in annual survival of PBs and recruit- 
ment were consistent (additive for annual survival) and 
of similar magnitude; about 0.25 higher for high than 
low-quality individuals at mean values, perhaps reflecting 
similar sources of variation. Lower survival and recruit- 
ment of young brant led to proportionally fewer low- 
quality individuals alive and reproducing in the older age 
classes and therefore limited our ability to identify 
remaining heterogeneity in older age classes. However, 
past and more directed studies of survival and recruit- 
ment indicated that condition at fledging affects future 
survival and recruitment (Sedinger and Chelgren 2007; 
Sedinger et al. 2008). We did not attempt to model 
sources of variation in clutch size of BRs; however, the 
range of variation we considered is consistent with past 
studies (Nicolai and Sedinger 2012). 

The effects of individual heterogeneity on population 
dynamics were tempered by how we parameterized heter- 
ogeneity in our models and also by our assumption about 
inheritance of quality. Although we did not observe heter- 
ogeneity in survival of adult brant, our results were sensi- 
tive to relatively small changes in the clutch size of adults 
because this class of individuals had the highest reproduc- 
tive values. However, differences in reproductive values 
among, and within, classes of individuals were highly 
dependent on assumptions about inheritance of quality 
class, and these differences were minimized or eliminated 
when inheritance was removed. Based on genetic studies 
of related species (e.g., Lesser Snow Geese, Chen caerules- 
cens caerulescens, Cooke et al. 1995) and our understand- 
ing of how brant goslings use brood-rearing areas 
(Lindberg and Sedinger 1998; Nicolai and Sedinger 2012), 



we suspect that quality is heritable either through direct 
genetic inheritance or through fidelity to high-quality 
environments or both. Therefore, harvest scenarios that 
included some level of inheritance of quality seem more 
appropriate than those that did not. 

Harvest dynamics of brant 

We demonstrated that individual heterogeneity alone 
could explain the observed relationship between harvest 
and population dynamics (Johnson et al. 1984; Lebreton 
2005). Under the scenarios we considered, compensation 
through heterogeneity resulted in about a 3% higher X at 
the highest rates of harvest compared with an additive 
scenario. Although current rates of fall harvest for brant 
have declined to about 1% (Sedinger et al. 2007), total 
harvest (fall and spring) may still have effects on popula- 
tion dynamics because (1) density-dependent reductions 
in recruitment and abundance from degrading habitat 
(Sedinger and Nicolai 2011) have likely produced a popu- 
lation composed of proportionally more high-quality 
brant, (2) spring harvest has likely increased in recent 
years and has definitely become a larger component of 
the overall harvest, (3) spring harvest occurs mostly near 
breeding areas and likely includes mostly high-quality 
individuals because migration timing and routes of pre- 
and nonbreeders do not expose those groups to as much 
spring harvest (Sedinger and Nicolai 2011), and (4) in 
spring, a large part of low-quality individuals have been 
eliminated, due to the heterogeneity in first-year survival, 
and the differences in reproductive values between low- 
and high-quality individuals are smaller, being only the 
result of the difference in recruitment pace, and possibly, 
clutch size. Therefore, the potential for heterogeneity to 
compensate for harvest through differential removal of 
low-quality individuals is likely lower in the spring than 
the fall (Ward et al. 1997; Francis 2000; but see Sedinger 
and Nicolai 2011). However, we do not currently have 
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adequate information on spring harvest rates or heteroge- 
neity in reproduction of brant to effectively model effects 
of spring harvest. 

Heterogeneity and harvest dynamics 

The concepts of individual heterogeneity and harvest are 
not new to hunting regulations. Populations are frequently 
managed for different quality individuals. For example, 
harvest regulations for large mammals commonly include 
restrictions on animals with different size antlers or horns, 
which may be linked to reproductive value of those indi- 
viduals (e.g., Garel et al. 2007). Species with some pheno- 
typic expression of quality (e.g., antler size) may more 
readily lend themselves to harvest regulations that account 
for individual heterogeneity. However, this type of pheno- 
typic expression is more common for males than females, 
and in particular for polygynous populations. Harvest of 
males has less effect on population dynamics than that of 
females. Nonetheless, excessive harvest of males with some 
traits can lead to undesirable population effects (Coltman 
et al. 2003). Gender-specific regulations may therefore be 
one way to minimize the effects of harvest in the presence 
of individual heterogeneity in quality. A more general strat- 
egy is to target different quality classes when they are tem- 
porally or spatially segregated, and this may be an 
appropriate strategy for migratory bird harvest. 

For migratory birds, most studies indicate that lower- 
quality individuals are more vulnerable to harvest (e.g., 
Greenwood et al. 1986; Dufour et al. 1993; Heitmeyer 
et al. 1993; Pace and Afton 1999) than are high-quality 
individual, but vulnerability is temporally and spatially 
variable. The potential compensating effects of heteroge- 
neity and differential harvest likely decrease during the 
season as both the proportion of low-quality individuals 
in the population and differences in reproductive values 
between quality classes decline. Spatially, we think that 
the most likely reason for variation in vulnerability is 
related to timing of migration for low- and high-quality 
individuals and is rooted in the spatial segregation of 
low- and high-quality individuals at the end of the breed- 
ing season because of different molting strategies/regions 
(Hohman et al. 1992). The timing and pattern of migra- 
tion may therefore differ because of variability in feather 
development, body condition, and local weather and hab- 
itat conditions. Low- and high-quality individuals may 
remain segregated during fall migration, and the extent 
and length of this segregation could vary annually. How- 
ever, harvest regulations are generally based on calendar 
dates and in any given year harvest may impact individu- 
als of different qualities, which may explain yearly varia- 
tion in the effects of harvest on populations. The 
compensating effects of heterogeneity may be maximized 



by considering both the timing and location of harvest 
(Francis 2000; Conroy et al. 2005). For waterfowl specifi- 
cally, we suggest that an individual heterogeneity model 
be explored as an alternative to the density-dependent 
models under consideration in current harvest manage- 
ment frameworks (e.g., Johnson et al. 1993). 
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